Descargate este archivo R-Markdown. Este archivo contiene todo el código de esta página.

Organizando tus librerias

Para este tutorial necesitaremos los siguientes paquetes. Importante: El paquete RevGadgets esta disponible directamente con install.packages(“RevGadgets”) pero vamos a revisar una versión especial que es la que grafica mapas estocásticos. Por esa razón vamos a bajar la versión que se indica aquí, así que. sigue las instrucciones como se indican abajo.

Convergencia del MCMC

  1. Baja los siguientes archivos y pegalos en la carpeta hisse_visualizaciones:

Revisemos la convergencia básica

# Con el paquete coda
mcmc_run1 <- readTrace(path ="hisse_polinizacion_run_1.log", burnin = 0.1)
## Reading in log file 1
mcmc_trace <- as.mcmc(mcmc_run1[[1]])
traceplot(mcmc_trace)

effectiveSize(mcmc_trace)
##              Iteration              Posterior             Likelihood 
##                 0.0000              1570.0326              4876.6295 
##                  Prior          extinction[1]          extinction[2] 
##              1551.4522               555.6097               439.0926 
##          extinction[3]          extinction[4]    extinction_alpha[1] 
##               665.4208               712.8727               484.9212 
##    extinction_alpha[2]     extinction_beta[1]     extinction_beta[2] 
##               305.5809               360.2756               234.7382 
## net_diversification[1] net_diversification[2] net_diversification[3] 
##              2265.9664               637.0799              6278.9435 
## net_diversification[4]                 q_0A0B                 q_0A1A 
##              5897.2707              7155.9250              4266.2468 
##                 q_0B0A                 q_0B1B                 q_1A0A 
##              3418.0186             18113.4734              3455.6445 
##                 q_1A1B                 q_1B0B                 q_1B1A 
##              3693.0056              5912.6036              3784.2706 
##    root_frequencies[1]    root_frequencies[2]    root_frequencies[3] 
##             24185.6143             73529.7285             72125.3665 
##    root_frequencies[4]          speciation[1]          speciation[2] 
##             71765.3892               804.1174              1036.4692 
##          speciation[3]          speciation[4]    speciation_alpha[1] 
##              3835.6708              4660.1254               853.2841 
##    speciation_alpha[2]     speciation_beta[1]     speciation_beta[2] 
##              1076.8263               877.0545               472.4078

Ahora hagamos un mejor trabajo con los gráficos de convergencia

# Esta parte se hace con ggplot2 y con RevGadgets
mcmc_run1 <- readTrace(path ="hisse_polinizacion_run_1.log", burnin = 0.1)
## Reading in log file 1
mcmc_run1<-data.frame(mcmc_run1)
mcmc_run1<- cbind(mcmc_run1,run=rep("run 1",length(mcmc_run1$Iteration)))

mcmc_run2 <- readTrace(path ="hisse_polinizacion_run_2.log", burnin = 0.1)
## Reading in log file 1
mcmc_run2<-data.frame(mcmc_run2)
mcmc_run2<- cbind(mcmc_run2,run=rep("run 2",length(mcmc_run2$Iteration)))

mcmc_table<-rbind(mcmc_run1,mcmc_run2)

trace_plot<- ggplot(mcmc_table, aes(x=Iteration,y=Posterior,group=run))+
              geom_line(aes(color=run))+
              theme_classic()
trace_plot

#Aqui se nota que hay que cortar un poco mas de burn-in (15,000 generaciones), asegúrense siempre de cortar mas valores si hace falta

Estimadores de las tasas

Queremos graficar y obtener estadísticas resumen de los parametros de interés. Recordemos que para el modelo hisse los parámetros que estimamos van a ser los siguientes:

Para los modelos más complicados siempre prefiero graficar violines, y tener más control en los colores y la figura, por eso mi selección de ggplot2. La comparación de violines es super útil.

Graficando tasas de transición

# Esta parte se hace con ggplot2
traitcols <- c("#3D348B", "#5A4FCF","#F18701", "#F35B04", "#1B7F3A", "#4CAF72","#D4A300", "#F2D15C")

hisse<- read.table("hisse_polinizacion_run_1.log", header=TRUE)
hisse<- hisse[-seq(1,15000,1),] # nunca olvidemos quitar el burn-in con esta instrucción tengo mas control de cuantas iteraciones remuevo.

tasas_transicion <- data.frame(dens=c(
         hisse$q_0A1A,hisse$q_1A0A, 
         hisse$q_0B1B,hisse$q_1B0B, 
         hisse$q_0A0B,hisse$q_0B0A,
         hisse$q_1A1B,hisse$q_1B1A),
  rate=rep(c("q_0A1A","q_1A0A","q_0B1B","q_1B0B","q_0A0B","q_0B0A","q_1A1B","q_1B1A"),each=length(hisse$q_0A1A))
  )

tasas_transicion$rate <- 
  factor(tasas_transicion$rate,
         levels = c("q_0A1A","q_1A0A","q_0B1B","q_1B0B","q_0A0B","q_0B0A","q_1A1B","q_1B1A"))


violin_transitions<- ggplot(tasas_transicion,aes(x=rate,y=dens, fill=rate))+
  geom_violin(trim=FALSE)+
  labs(title="Tasas de transicion")+
  scale_fill_manual( values = traitcols)+
  theme_classic()
violin_transitions

¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta

Ahora vamos comparar Transiciones entre estados escondidos

## Vamos a crear una tabla de las diferencias entre las transiciones
test_tasas_transicion_hisse<-
  data.frame(
  dens=c((hisse$q_0A1A - hisse$q_1A0A),(hisse$q_0B1B - hisse$q_1B0B),
         (hisse$q_0A1A - hisse$q_0B1B),(hisse$q_1A0A - hisse$q_1B0B),
         (hisse$q_0A0B - hisse$q_0B0A),(hisse$q_1A1B - hisse$q_1B1A),
         (hisse$q_0A0B - hisse$q_1A1B),(hisse$q_0B0A - hisse$q_1B1A)),
  difference=rep(c("q_0A1A - q_1A0A","q_0B1B - q_1B0B",
                   "q_0A1A - q_0B1B","q_1A0A - q_1B0B",
                   "q_0A0B - q_0B0A","q_1A1B - q_1B1A",
                   "q_0A0B - q_1A1B","q_0B0A - q_1B1A"
                   ),each=length(hisse$q_0A1A)))

test_tasas_transicion_hisse$difference <- 
  factor(test_tasas_transicion_hisse$difference,
         levels = c("q_0A1A - q_1A0A","q_0B1B - q_1B0B",
                    "q_0A1A - q_0B1B","q_1A0A - q_1B0B",
                    "q_0A0B - q_0B0A","q_1A1B - q_1B1A",
                    "q_0A0B - q_1A1B","q_0B0A - q_1B1A"))

violin_hidden <- ggplot(test_tasas_transicion_hisse,aes(x=difference,y=dens, fill=difference))+
  geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+ #Intervalos de maxima probailidad
  labs(title="Tasas de transición entre estados escondidos ")+
  geom_hline(yintercept = 0, linetype = "dashed", lwd = 1)+
  scale_fill_manual( values = traitcols)+ 
  theme_classic()+
  theme(legend.position = "none",  plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 45, hjust = 1))
## Warning: The `draw_quantiles` argument of `geom_violin()` is deprecated as of ggplot2
## 4.0.0.
## ℹ Please use the `quantiles.linetype` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
violin_hidden

¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta

Graficando tasas de diversificación

# Esta parte se hace con ggplot2

divcols<-c("#E63946","#1D3557","#F3A5AB", "#457B9D")

# Recuerda de RevBayes 1=0A, 2=1A, 3=0B, and 4=1B
tasas_diversificacion_neta <- data.frame(dens=c(hisse$speciation.1.-hisse$extinction.1.,hisse$speciation.2.-hisse$extinction.2., hisse$speciation.3.-hisse$extinction.3.,hisse$speciation.4.-hisse$extinction.4.) ,rate=rep(c("net_div_0A","net_div_1A","net_div_0B","net_div_1B"),each=length(hisse$speciation.1.)))

tasas_diversificacion_neta$rate<- factor(tasas_diversificacion_neta$rate,levels=c("net_div_0A","net_div_1A","net_div_0B","net_div_1B"))

violin_diversificacion<- ggplot(tasas_diversificacion_neta,aes(x=rate,y=dens, fill=rate))+
  geom_violin(trim=FALSE)+
  labs(title="Tasas de diversificación neta")+
  scale_fill_manual( values = divcols)+
  theme_classic()
violin_diversificacion

¿Qué podemos concluir de este gráfico?

¿Prueba de hipótesis: La polinización esta asociada a la diversificación?

Una de las ventajas más importante de las distribuciones posteriores es que rápidamente podemos concluir si dos parámetros son iguales o distintos. Formalicemos entonces la prueba de hipotésis para saber si el caracter esta asociado a la diversificación.

Construimos dos estadísticas de resumen que son las diferencias entre la diversificación neta de 0 y 1 para el estado A y el estado B respectivamente:

Si estas diferencias son 0 con alta probabilidad esto significa que la polinización y sus estados no hacen la diferencia en el proceso de diversificación.

# Esta parte se hace con ggplot2

difcols<-c("#FF006E","#FFC2DC")
T_diferencias<- data.frame(dens=c((hisse$speciation.1.-hisse$extinction.1.)-(hisse$speciation.2.-hisse$extinction.2.),(hisse$speciation.3.-hisse$extinction.3.)-(hisse$speciation.4.-hisse$extinction.4.)),diferencia=rep(c("T_A","T_B"),each=length(hisse$speciation.1.)))

violin_difference<- ggplot(T_diferencias,aes(x=diferencia,y=dens, fill=diferencia))+
  geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
  labs(title="Estadisticas resumen")+
  scale_fill_manual( values = difcols)+
  geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
  theme_classic()

# Observemos que 0 cruza la diferencia entre las tasas lo que quiere decir que pueden ser iguales
violin_difference

¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta

Formalizando la prueba de hipótesis \[ H_0: T_A=0 \textrm{ y } T_B=0\]

# Esta parte se hace con dplyr

cuantiles <- T_diferencias %>% group_by(diferencia)%>%reframe(res=quantile(dens,probs=c(0.025,0.975)))
cuantiles
## # A tibble: 4 × 2
##   diferencia     res
##   <chr>        <dbl>
## 1 T_A        -0.284 
## 2 T_A         0.390 
## 3 T_B        -0.636 
## 4 T_B        -0.0135

Observamos que el intervalo de credibilidad del 95% de la distribución de \[T_A\] es (-0.284, 0.389) y el intervalo de credibilidad del 95% de la distribución de \[T_B\] es (-0.636, -0.013). Cómo 0 pertenece al intervalo \[T_A\], significa que el tempo de la diversificación para A es igual, mientras 0 no pertenece al intervalo \[T_B\]entonces significa que el tempo de la diversificación para el estado 1 es diferente. Importante: Noten que en las pruebas de hipótesis bayesiana no utilizo, p-valores, ni significancia, ni rechazo, solo probabilidad e intervalos de credibilidad. Esto sucede porque la manera en que interpretamos probabilidad en bayesiana es distinto. Tengan cuidado cuando interpreten.

Segunda parte: ¿Qué pasa con los estados escondidos?

Construimos dos estadísticas de resumen que son las diferencias entre la diversificación neta de A y B para el estado 0 y el estado 1 respectivamente:

Si estas diferencias son 0 con alta probabilidad esto significa que hay cambios en el tempo de la diversificación pero se deben a algo más que nosotros no medimos en la filogenia.

# Esta parte se hace con ggplot2

difcols<-c("#1DB32C","#BFEEC3")
T_diferencias<- data.frame(dens=c((hisse$speciation.1.-hisse$extinction.1.)-(hisse$speciation.3.-hisse$extinction.3.),(hisse$speciation.2.-hisse$extinction.2.)-(hisse$speciation.4.-hisse$extinction.4.)),diferencia=rep(c("T_0","T_1"),each=length(hisse$speciation.1.)))

violin_difference<- ggplot(T_diferencias,aes(x=diferencia,y=dens, fill=diferencia))+
  geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
  labs(title="Estadisticas resumen")+
  scale_fill_manual( values = difcols)+
  geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
  theme_classic()

violin_difference

¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta

Formalizando la prueba de hipótesis \[ H_0: T_0=0 \textrm{ y } T_1=0\]

# Esta parte se hace con dplyr

cuantiles <- T_diferencias %>% group_by(diferencia)%>%reframe(res=quantile(dens,probs=c(0.025,0.975)))
cuantiles
## # A tibble: 4 × 2
##   diferencia    res
##   <chr>       <dbl>
## 1 T_0        -0.453
## 2 T_0        -0.195
## 3 T_1        -1.02 
## 4 T_1        -0.339

Observamos que el intervalo de credibilidad del 95% de la distribución de \[T_0\] es (-0.452, -0.195) y el intervalo de credibilidad del 95% de la distribución de \[T_1\] es (-1.024, -0.339). Cómo 0 no pertenece a estos intervalo entonces \[P(H_0\lvert Datos)<0.05\] lo que significa que el tempo de la diversificación para A y B es diferente. Este resultado unido al resultado anterior indican que el tempo cambia debido a algo más que nosotros no medimos pero se manifiesta en el árbol. Estos resultados son equivalentes ajustar el modelo de caracter independiente CID y seleccionarlo, indicando que la diversificación es independiente del caracter de interés.

¿Qué hubieras concluido con BiSSE?

Por que la curiosidad mató al gato chequemos que paso con BiSSE

Baja los siguientes archivos:

# Esta parte se hace con ggplot2
bisse<- read.table("bisse_polinizacion_run_1.log", header=TRUE)
bisse<- bisse[-seq(1,15000,1),] 
divcols<-c("#E63946","#1D3557")

# Recuerda de RevBayes 1=0A, 2=1A, 3=0B, and 4=1B
netdiversification_rates<- data.frame(dens=c(bisse$net_diversification.1., bisse$net_diversification.2.) ,rate=rep(c("net_div_0","net_div_1"),each=length(bisse$net_diversification.1.)))

violin_diversification<- ggplot(netdiversification_rates,aes(x=rate,y=dens, fill=rate))+
  geom_violin(trim=FALSE)+
  labs(title="Tasas de diversificación neta")+
  scale_fill_manual( values = divcols)+
  theme_classic()
violin_diversification

Prueba estadística

difcols<-c("#1DB32C","#BFEEC3")
T_diferencia<- data.frame(dens=(bisse$net_diversification.1.-bisse$net_diversification.2.),diferencia=rep("T",each=length(bisse$net_diversification.1.)))

violin_difference<- ggplot(T_diferencia,aes(x=diferencia,y=dens, fill=diferencia))+
  geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
  labs(title="Estadìstica resumen")+
  scale_fill_manual( values = difcols)+
  geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
  theme_classic()

# Observemos que 0 cruza en una probabilidad muy pequeña de la distribución la diferencia entre las tasas lo que quiere decir que pueden ser iguales
violin_difference

## 0 no esta en el intervalo de probabilidad del 95%
quantile ((bisse$net_diversification.1.-bisse$net_diversification.2.),probs=c(0.025,0.975))
##      2.5%     97.5% 
## -0.558164 -0.092405

¿Qué hubieras concluido para BiSSE? ¿Cuál es la mayor falla de este modelo?